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Introduction 


In a recent paper Phillips and Rose [2] described a compact system of 
finite difference equations to treat the equilibrium of elastic bodies. For 
any brick cell of volume 0(h 3 ) within the body these equations express, to 
second order accuracy, the relationship between the traction forces and the 
displacements on the faces of the cell which results when the cell is in 
isolated equilibrium. Global equilibrium then occurs when the net traction 
force across a face common to any two cells vanishes; this condition serves to 
eliminate the traction forces as variables and results in algebraic conditions 
for equilibrium which are expressed in terms of the displacement values on the 
sides of neighboring cells. These algebraic equations seem particularly 
suited to solution by iterative techniques which can exploit the parallelism 
of computer architectures. 

In this method, since each cell is in isolated equilibrium, the potential 
energy in the cell is equal to half the work done by the traction forces and 
displacements on the surface of the cell. The total potential energy within 
the body can then be calculated, and the principle of minimum potential energy 
simply results in the condition that the net traction forces vanish across 
common cell faces. 

This paper examines this idea when general volume cells are employed. An 
approximate equilibrium condition for each cell is obtained by means of the 
construction of a transmission matrix which, when applied to displacement 
values on the faces of the cell, yields approximate traction forces on cell 
faces thus producing equilibrium in the cell. We call these equations compact 
finite element equations . Equilibrium throughout the body is then obtained by 
determining displacement values on the faces of cells such that the net 
traction force vanishes on every face incident to two cells. 
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The transmission matrix arises as an approximation to a boundary integral 
operator on a cell. As we shall see, it is also closely related to the 
ordinary stiffness matrix in finite element methods. 

We are content here to describe a method of construction which yields 

2 

0(h ) convergence. This is done in Part I, which shows the equivalence to a 
special nonconforming finite element method. The method can easily be 
converted to a stress, rather than a displacement, formulation. 

Part II discusses results which have a bearing upon the practical 
implementation of the general construction. Using nonconforming linear 
elements the compact equations for a tetrahedron are easily obtained. We show 
that the compactness property is preserved under a condensation technique; as 
a result, general cells can be handled without constructing an enlarged 
approximation basis in such cells which the general theory would require. 

In the context of the results of [2] the practical implication of these 
results in three dimensions is: in order to approximate the displacements of 
a body occupying a volume ft, we need only cover ft by regular brick cells 
and introduce tetrahedral cells near the boundary of ft. The compact 
difference equations for brick cells and the compact finite element equations 
for tetrahedral cells then establish the equilibrium conditions for individual 
cells in ft. Global equilibrium in ft is then obtained by imposing the 
algebraic balance of traction conditions throughout ft. The general problem 

of developing effective iterative methods for this type of construction will 
be examined elsewhere. 

The paper concludes by suggesting an approach to the homogenization 
problem using the condensation technique mentioned above. 


Part I 


Equilibrium of Linear Elastic Bodies 

Notations: 

In the following x = (x^x^Xj) is a point in I? and fl is a bounded 
volume with boundary T. {w} M indicates a partition of fl into M 
subvolumes, or cells; y = Y(<*>) is the boundary of to on which n is a unit 
outward normal, y = y(to,to') indicates the surface incident to both to 
and to' and is called an interior face; faces of y(to) which lie in T are 

a 

denoted by Yp* 

Aw denotes the volume of 0) and Ay the surface area of Y* We assume 
positive constants cg^ which are independent of M such that if 
h = 0( 1/M) 

0 < Cq h^ _< Aco £ c^ h^. 

Finally, if v is a continuous function on fl then v denotes the 
values of v on a cell boundary y. 

1.1. General Results 

In the following B is a material body occupying a volume fl; 
u = displacement, t = stress tensor, and e = strain tensor. Linear theory 
leads to 

e(u) = grad(u^ + u C ), (1.1a) 


t(u) = t(e(u)); 


(1.1b) 
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in which (b) represents the constitutive relationship of the material. Both 
e and t are symmetric, and the reciprocal principle of Betti-Rayleigh 
(Timoshenko and Goodier [4]) is expressed by 

W(u,w;fl) = V 2 I e(u)x(w)du = W(w,u;Jl). 

Q ~ 

Let 

_ o 

Lu = div t(£) - q u, 

where q is a nonzero constant; the case q = 0 will be treated later. Then 
the boundary value problem 

Lu = 0 in Cl, 

" = iir on F (1.4) 

describes the equilibrium of the body subject to a spring load. The potential 
energy is given by 

U(ti;fl) = W(u; Cl) + V 2 ^ dto. (1.5) 

n 

The principle U(v;J2) = min then leads to the solution u of (1.4) and forms 
the basis for many finite element methods. 

The solution of (1.4) is determined solely by the boundary data and may be 
represented by a solution operator E as 

u(x) = E(x)up 


( 1 . 2 ) 


(1.3) 


( 1 . 6 ) 


in which E(x) reduces to the identity operator for xeT. The traction force 


on T resulting from a displacement 11 is given as 

p = t(u) * n| p. 


(1.7) 


For the equilibrium solution (1.6) we may write 

p = p(u_p) = x(Eup) *ri| p = Tp u_pj 


( 1 . 8 ) 


we call the boundary operator Tp a transmission operator . 

For arbitrary displacements u,u, standard integral theorems are 

^ u^ pdy = 2U(u;fi) + / u C Lu dio 

r — 1 ft 


and, using (1.2), 


£ (uH p - Up p)dy = / (u* 1 Lu - u C Lu)du>. 

r ft 


At equilibrium, Lu. — 0; thus 


$ Up £(u r )dY = 2U(u;ft) 


(1.9) 


f p(u.iO dY = 4 nr £(nr^ dY * 
n r 


and 


( 1 . 10 ) 
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Equation (1.9) simply states that at equilibrium the work due to the traction 
forces and displacements on T equals twice the potential energy in fi; 
equation (1.10) expresses the reciprocal principle (1.2) comparing two 
equilibrium states. 

A representation formula for the transmission operator Tp can be 
obtained by setting 

a 

u.rO^') = <5(x - x') , x,x' e r (1.11) 

in (1.10) where 6 is the Dirac function. Defining 

T r (x,x') = £(u.rQi>x')) (1.12) 

(1.10) leads to 

p(x;u r ) = £ uf(x')T r (x,x')dy(x') = Tp(x)-(up). (1.13) 

It follows that Tp is a symmetric operator. Also, since 

2U(_u;fl) q 2 / _u 2 dm > 0 (1.14) 

Q 


(1.9) results in 


so that Tp 


lj> 1L J>(u.iO dY = 4 Up Tp(up)dy > 

is positive definite when 


q 2 / u 2 do» > 0 (1.15) 

> 0. We may summarize these 


results in 
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Theorem 1: Let 8 be in equilibrium due to prescribed surface 

displacements Up. Then the resulting traction forces on T are determined 
from u r by means of a symmetric positive definite transmission operator Tp 
by ja( Up ) = Tp(up), and the work done by these traction forces is equal to 
twice the potential energy, i.e. , 

<f> Up p(u r )dY = 4 Hr T r(— = 2U(u;£l) (1.16) 

pi r 


where _u = Eup . 

1.2. The Equilibrium Problem with Constraints 

Introduce a partition of £1 onto M volume cells {w}^. Let 1/ denote 
the class of continuous functions v on £1 which satisfy the boundary 

conditions v = u p on T. The values of v on the boundary surfaces y(o)) 
of cells 0 ) are v^^ or, simply, v y * 

For a given vel/ consider the problem of minimizing the potential energy 
of the system by displacements v/ subject to the constraints Wy(oj) = — y((o)’ 
for aie£2, i.e., 

min U(w;£l) > U(Eu r ;£l). (1.17) 

w 

The solution is given by w = Ev.^^ , ue£l, and in each cell 

min 2U(w;u>) = 2 u(Ev y(w) = cf vj pCv^dY. (1.18) 

- Y(w) 

Denote the total work due to traction forces against the constraints by 
Q(v;£l). Then, using (1.13), 


- 8 - 


Q(v;«) = l v y P(v )dy = l v* T (v )dy > 2U(u;B) >0. (1.19) 

weft ' ' aiefj t t — y — — — 

Any two displacement constraints v>v'eV can be considered equivalent if 

-Y(oj) = — y((o) ’ weft > since they lead to identically equal equilibrium 
displacements in each cell, i.e., w - w' = E(v y - v') = 0. The consequences 
may be stated as 

Theorem 2: The total work Q(_v; f2) due to traction forces 

pCv^) = T y (£ y ) and the constraints v y is minimized by the class of 
constraints vet/ which satisfy v y(u) = u y(aj) where u_ is the equilibrium 
displacement of problem (1.4). 

The net traction force acting at a point 1 on a face y common to 
neighboring cells w,(o will be denoted by 

l.(i>v) = p(i,v Y ) + £ (i,v r ) = T y (i)v y + T r (i)v r . (1.20) 

Also we will let £ r = £ r (v y ^) indicate the traction force on the face of a 
cell to which lies on T, Then summation by parts yields 

q -• a> ' 2> tm * £< -^ ,dT -f ****** 4 { ^ ~y d *• (1 - 21> 

Obviously, the conditions 

f *d*v) = 0, yeft (1.22) 

Y 

are the Euler conditions for Q(v;ft) = min. The result is 
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Theorem 3: A necessary and sufficient condition that Vy( M ) = » 

weft, is that the net work done by traction forces vanishes on every interior 
face on which constraints are imposed; at equilibrium, £p = p(£p) on ^ 
( 1 . 21 ). 

Introduce the norm 

II v I'll = (b I / v^vdy) ^2 . (1.23) 

A A 

yeft y 

If w = Ev^ , 

/ v C T (v )dy = 2 U(w;oj) q / w wdio 
y(w) Y Y -Y a) 

= aq 2 h j v y V v dY 

Y(u) _Y Y 

O 

where a = 1 + 0(h ). Thus a lower bound for the lowest eigenvalue of 

the transmission operator T^, is 

A = ahq 2 . (1.24) 

An estimate for the difference between the constrained and unconstrained 
problems can now be obtained: let _e = _v - u^ so that Cp = 0. Then 

(1.21) yields 


Q(£,ft) = l $ T Y (£ )dy = l dy (1.25) 

weft y(w) a a y y 

y eft y 

in which _f a = f . (e^ , e^ ) so that 
Y Y 
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i.e. 


y 


aq 2 leljj £ h_ 1 le*h "f(£)« h , 
oq 2 lel h < h _1 Ilf (e) n h . 


1.3. An Approximation Method 

In this section we describe a method for approximating a smooth solution 
<£ of (1.4) with O(h^) accuracy. 

Consider a typical cell u whose boundary y(w) consists of k faces 

A A A 

Y 1 ’^2’*** » Y k* Reca11 the representation (1.13) of the transmission operator 
Ty as the integral operator 

T Y QL;£ Y > = $ T Y (x,QvU)dY(£) 

in which T (x,£) is the traction force due to a displacement 

5(|x. “ £|) on Y • Let £ indicate the centroid of a face y whose area is 

A 

Ay . Define 

[Z-yT H [v(£i )Ay 2 ,v(i 2 )AY 2 , * * ‘ , y ( J » (1.27a) 

k 

T*?(x)[v y ] = l T Y (x,i i )v(£ i )AY i> (1.27b) 

[TyJ = (X Y Cl 1 5i. j >)- (1.27c) 

The kxk matrix [t ] is symmetric and positive definite and 


£. h (i>2L Y ) = T^}(i)[v Y j 


(1.28) 
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approximates the traction force on a side Y* 

Q h (v;“) - ][*,<„>! (1 - 29) 

approximates the work done by traction forces and displacements on y(o>). The 
net traction force on a face y common to adjacent cells w and u> is then 
approximated by 

fj(5 , [vj ) = T^(i)[v Y j + T^(i)[v r J = £ h (i,v y ) +£ h (i,v Y .). (1.30) 

Y 

Let us agree to approximate the boundary condition v^, = iip for vel/ by 

Ay r »v(£ r ) = / u p dYp, Yp e r (1.31) 

Y r 

in which tip indicates the centerpoint of a face Yp of a cell which lies on 

r. 

We propose to show that the solution u* 1 of the following algebraic 
problem P* 1 provides an O(h^) approximation to the solution u. of (1.4) at 

A A 

points 5 of Y e ^« 

Problem P h : Solve 


f^(£>[v]) - 0, Y^n 

Y 



subject to the boundary conditions (1.31). 


(1.32) 


for 
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First, using (1.29), define 


Qk(v;&) = J Q^Cvjid). 
ueJl 


Then summation by parts leads to 


(1.33) 


Q h (v;Sl) = l [y y ] t [Tjj][v ] = l v t (5 r )p h (L)A Y 
meft ITT i — 1 T 

C r er 

+ l V A (i)fh(£,[v])Ay (1.34) 

* Y Y 

yen 

in which p (Ej.) represents an approximate traction force arising from a cell 
which has a face lying on T. 

The value V 2 Q (v;ft) given by (1.33) is an approximation to the potential 
energy of the constrained system. The problem of minimizing Q h (v; SI) among 
the class of algebraic constraints satisfying (1.31) leads, using (1.34), to 
(1.32) as Euler conditions. 

Next, introduce the norm 


5 (h l vHbzl i)A?(|))V2. (1.35) 

left 


From earlier results 


Q h (v;Sl) > aq 2 [[v]£ 


(1.36) 


where a = (1 + 0(h 2 )). 

With u indicating the solution of (1.4), define 
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e = u - v. 


(1.37) 


Then e satisfies the boundary conditions e(5_ r ) = 0* Also, from (1.32), 

l h (i;[e.]) =ll*(i;[u]). Hence (1.34) and (1.36) lead to the estimate 
Y Y 

<*q 2 [| ejl h < h” 1 [| 1*Ir] |] h . (1.38) 

Y 

Using (1.27) and (1.13), define 

(p h - p)(u Y ) = T^j = T Y (u y ). (1.39) 

Then, since f_~ (u.) = 0 
Y 

fj[u] = [f* - f J(u) = (p h - P)(u y + HyO (1 *40) 

Y Y Y 

where y is the face common to y and f ' • 

2 2 

For illustration, consider a boundary value problem for V u - q u = 0. 
In this case both the solution u and the kernel of the transmission operator 
can be represented as a superposition of exponential solutions of the form 
exp q(£*x). Then 

|(p h " P)(u Y )| = q 2 0(h 3 ). 


A similar superposition principle applies to the problem (1.4) when, as we may 
assume, the material properties do not vary in a cell. Using (1.40) we then 


obtain the estimate 
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[li^u) |Jh = q 2 o(h 3 ) (i.4i) 

y 


and (1.38) yields 


Thus 


[|e|] h = 0(h 2 ). 


(1.42) 


Convergence of P : The solution u^ 1 of problem P^ 1 converges to the 

solution u of (1.4) with an accuracy 0(h 2 ) when u is sufficiently 
smooth. By continuity this remains true when q = 0. 


We may summarize matters as follows: for each cell which is in 

equilibrium assume that an approximate algebraic relationship can be 
established between the traction forces [pjjj and displacements [v y J on the 


faces of the cell by means of a transmission matrix 



viz. , 


[PyJ = [ T IjJ[v y J. (1.43) 

We call (1.43) a compact finite volume scheme . Next, establish global 

equilibrium in ft by imposing the balance of traction conditions 
ll 

JL„(v) = 0; this corresponds, at the algebraic level, to the elimination of 
Y 

stresses and strains in (1.4) so as to obtain a second order system of 
equilibrium of partial differential equations for the displacements. When the 
boundary conditions (1.31) are imposed, the resulting algebraic equations have 
a unique solution u h and [|u - u 1 *^ = 0(h 2 ) when u is smooth. 



-15- 


1.4. An Abstract Formulation 

2 

Consider again the equilibrium equations Lu = div t(u) - q u = 0. We 
may, in principle, determine a family of solutions i = l>2,***,k in 

any cell which has k faces such that 


Then 


Wit* 


ir 


k 

w(x) = l 0. (x,C. )v(£, ) 
i=l 


(1.44) 


provides a solution such that w(j^) = jv(j^) for JLi G Y±* I n an Y cell we 
may write 

w(x) = E(x)P h [v y J (1.45) 


where E(x) is the solution operator for the boundary value problem (1.4) 
and P* 1 is a projection operator. The traction forces which result are 

P h QL>[z. Y J “ T CEP h [v Y J ) -n| x = T^(x)[v y j. (1.46) 

Thus 

ll!Jj - (T^J (1.47) 


with r T^ 1 1 = P T^ 1 yields the required compact equations on the cell. 

L Y J Y Y 

This abstract formulation indicates how higher order approximations can, 
at least in principle, be obtained. We shall not pursue this matter here, 


however. 
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The construction (1.44) was proposed in [3] as a general method for 
obtaining the solution of elliptic partial differential equations. The 
present paper clarifies the physical ideas involved and presents the estimates 
required to establish a general convergence proof. 

1.5. Further Considerations 

The previous discussion was focused upon the Dirichlet problem for (1.4) 
in which displacement values Up were prescribed on T. The numerical 
treatment of this type of boundary condition is described by (1.31). Problems 
in which the traction force is prescribed on a part of the boundary r° , say 

jKO - p (£) > £ G F , can be treated by using (1.28) and imposing the traction 
condition in the form 

T^ 0 (l)[v y J = £ °(£), (1.48) 

A 

where £ indicates the centroid of the face of a cell boundary y° which 
lies on T®. 

Finally, instead of formulating the algebraic equilibrium equations in 
terms of the displacements by means of the balance of traction conditions 
— » ^ = t ^ ie tractions may be considered as the primary variables and the 

y 

problem solved in terms of them as follows: First note that the tractions 

[ By J as given by (1.43) are necessarily compatible with the displacements 
[■By -I* From 0*24), when q * 0, the minimum eigenvalue of the transmission 
matrix [T y ] was shown to approximate ahq 2 > 0. We may thus set 

[Syj = [t ijr 1 


and write 


(1.49) 
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* [ s ?l[£ Y i- 


a. 5o> 


Consider 


R h (p) - [£ Y ( 0) )J t [ S Y( u ,) J[£ Y (a))J* 


tuefl 


(1.51) 


Let us agree to use (1.50) to express a boundary condition of the form (1.31) 
in terms of [p^J. Consider the variational problem 


max R (p) 
— Y 


(1.52) 


where p satisfies the balance of traction conditions f* 1 = 0 across any 

£.y 

Y 

interior face Y and also satisfies boundary conditions on T in the form 
just indicated. A summation by parts procedure similar to that indicated by 
(1.34) now leads to Euler equations which express conditions that the jump in 

L A 

displacement values of v across any interior face y must vanish. Using 
(1.50), these conditions are expressible in terms of the traction variables on 
the faces of neighboring cells. The result (c.f. [4]) is simply a restatement 
of Castigliano's principle in discrete form. 


Part II 

Practical Developments 

We now turn our attention to placing the previously described theoretical 
development into a practical framework. Our preliminary discussion will 
concern the Laplace equation in two— dimensions; this will enable us to 
illustrate several key features in a more familiar context. A later section 
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describes the construction details required to treat the equilibrium problem 
for the elastic body described in Part I. 

Suppose is a domain in ft can be conveniently covered by 

rectangular cells as shown in Figure 1 . 



Figure 1. A covering of a domain ft by rectangular cells. 

Interior cells can be treated by compact schemes for rectangles; cells at the 
boundary can be reduced to the treatment of compact schemes for triangles. 

Compact schemes for rectangles have been developed using the finite 
difference calculus in a manner described by Phillips and Rose [1], [2]. We 
shall indicate below how these can be obtained by a simple condensation 
process applied to more elementary results for triangular cells which we first 


discuss. 
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II. 1. A Model Problem 

2 

In this section we treat the Laplace equation 7 u = 0 as 
problem. 

a) Triangular Cells 

Consider a triangular cell to whose area is Aw. Choose 

coordinate system having the centroid of the vertices , i= 
origin 0. The centerpoints of the sides are » 5.2 >1L3 > as s ^ own 


x 2 



5 


Figure 2. A triangular cell. 

Let £(x,5^), i=l,2,3, denote area coordinates for the triangle 

3 3 

x= I *(x,5.K., I *(x,5.) = 1. 

“ i=l ” i=l 

The vertices [£_] and [x.] are thus related by 


a model 

a local 
1,2,3, as 
Figure 2. 


( 2 . 1 ) 
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m = n[x], 



i 

i 

0 




( 2 . 2 ) 


Also, n(C^) = n., 1-1, 2, 3, indicates the outward normal vectors on the 
sides. Finally, Ay^ = “ 2ifJ I s the length of a side containing 

and 

mC^) = Ay i «n(^ i )/Aaj. (2.3) 

Noting our choice of origin, A(x,C) then has the form 


A (£»i.i ) = (j + x*m(£i)) • 


(2.4) 


The function 


3 

w(x) = l %I-)v(?.) 
i=l 1 1 


(2.5) 


2 

is a solution of V w = 0 for which w(^) = vU^). The 'traction' defined 
by 

P h (Ii>v) - £(!!>. Vw^) (2.6) 


is thus given by 
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p h (£ i ,v) = l t^C^JnC^KAu) 


so that 


[p?j - [t;j[v y j 


(2.7) 


is a compact scheme on w with the transmission matrix 


[t^J 5 l . 


( 2 . 8 ) 


Next consider two adjacent triangular cells u),a)^ having the side y in 
common (Figure 3). 


e 


3 




Figure 3. Two neighboring triangular cells a ),w'. 


The balance of traction condition at 5^ is 


f'dpv) = Pydpv) + p^dj.v) = o 


(2.9) 
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which reduces, using (2.8) to 


VaS + S=- )v( ii > = ‘ AS J 2 sf Si 

- A^ j <£i>‘ Si vCipAYl 


Suppose (D,co" are acute triangles; then 


Define 


—1 — i ^ i-2,3. 


“u s - Si Sj AYj/AY,, 


p = AO)'/ ( Ad) + AoT) , 


p" = Aw/ (Aw + AuT) . 


Then v^_. > 0 and (2, 10) leads to 


v< Il ) ■ J 2 (pU li v( ii ) + p ' v Ij v( £p) 


However, 


<j> ndy - 0 
Y(t») 


so that 


Thus 


v 12 + v i3 " - n5lH. 2 a ^ 2 + n 3 by 3 )/&y l = 1. 


I ,(.P v ii + P' v l ± ) = P + P' = 


i=2 


( 2 . 10 ) 

( 2 . 11 ) 

( 2 . 12 ) 

(2.13) 

(2.14) 

(2.15) 


1 


(2.16) 
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Hence v(5^) as given by (2.13) is an average with positive weights of the 
neighboring values of v in the triangles. 

Figure 4 illustrates two situations. 


2 3 ' 




Figure 4, Cases of two neighboring right triangular cells. 

In the situation in Figure 4a, (2.13) leads to 

vUj) »V4 (v(C 2 ) + vUp + v (jLp + v(5p) (2.17) 

2 

which is a familiar 5-point approximation to 7 u = 0. In the situation in 

Figure 4b, the coefficients of the terms v(j^)> v(Ep in (2.13) vanish so 


that 


v(£ x ) = V2 (v(|_ 2 ) + v (Ip) 


(2.18) 
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Hcre we observe a difference between the result in this paper and that of 
the typical finite element treatment of V 2 u = 0 (Zienkiewicz [5]). In the 

r\ 

latter a linear approximation to V u = 0 in u is conventionally 
described in terms of area functions Ux.x^ for the triangle using the 
vertex points {x^, viz., 


3 

w(x) = l L(x,x.)v(x ). (2.19) 

i=l 1 1 

This construction leads to conforming approximations, i.e., approximations 
which are continuous across the sides of neighboring triangles. In contrast, 
the approximation (2.5) is discontinuous across cell sides except at the 
center points (S^)> i.e., (2.5) is nonconforming. For conforming 

approximations the standard 5-point finite difference formula expresses the 
relationship between a vertex value and its four neighboring vertex values, 
shown in Figure 5, viz., 

V 1 =1/ * (v 2 + V + v 3 + v 3' ) (2.20) 



2 ' 


Figure 5. The relationship of vertices involved in the 5-point 
formula (2.20) using conforming approximations in cells. 
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This example is often used to point to the close relationship between 

finite element and finite difference formulations of this problem. However, 

the theory given in Part I shows that the scheme represented by (2.17) and 

(2.18) also converges with O(h^) accuracy and thus represents a valid, if 

2 

less conventional, finite difference approximation to V u = 0. 

The next example will help illustrate an important property of 
approximations based on the nonconforming approximation (2.5) which is not 
shared by the conforming approximation (2.19). 
b) Rectangular Cells 

2 

We propose to show that, for V u = 0, compact schemes for rectangular 

cells can be developed by one of two equivalent means. One uses the general 

approximation method described 'in Part I. The other uses a condensation 

process based upon the example just given for triangular cells. Both lead to 

the same result. For reasons to be described, the condensation technique 

fails when the conforming approximation (2.19) is employed. 

i) An application of the previous theory: 

2 2 2 

The functions (l,x,y,x - y ) are solutions of V u = 0. From these we 
may construct solutions such that 

3 

w(x) = l 8(x;C,)v(?. ) (2.21) 

i=l 1 1 

2 

is a solution of V u = 0 in a) such that w(j^) » v^^)* 

2 

Write V u = 0 in system form as 

r + s =0 
x y 


r 


u , 
x* 


y 


s = u 


( 2 . 22 ) 
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Figure 6 illustrates a square cell on whose sides the indices 
(i ± l,j ± 1) indicate points corresponding to the centerpoints of sides of 
length 2h. 


<U+l> 



0+1, j) 


Figure 6. A square cell having indical labeled sides of length 
2h based upon a triangular decomposition. 

The 'traction' values at the centerpoints of the sides, labeled 
counterclockwise, are 


3 K,j+l> “ r i-l , j ’ " 8 i,j-l» r i+l,jJ ' 


( 2 . 23 ) 


Vi h 

The transmission matrix [t^J provides the relationship between [p^J and 


[« y J* 


Introduce the standard finite difference notations 
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5 

x 


i* j 


(r 


i+l,j 


r 


i-l.j 


)/2h, 


w x r i,j <r i+l,j + r i-l,j )/2 

with corresponding definitions for 5 s. .,y s. ,. In [1] the following 

y 1 > j y 

compact scheme was described: 


6 r . . + 6 s. . = 0 
x i,j y i»J 


“x r i,j ■ 5 X “y s i,J ' S y "i,r 


,2 h 2 

h r h . 

u r. . 7 t o r. . = p s. . 7? o s. . 

x i,j 2 x i, j y i,j 2 y i,j 


(2.24) 


The fact that the approximation (2.21) leads to (2.24) may be verified by 

2 2 

noting that each of the functions [l,x,y,x - y j in the approximation basis 
satisfies these equations. We may write (2.24) as 

r i±l,j ■ K 4h ' 1( “x - V Ju i,j 

s i,;i±l - [ { y 4 h ’ 1( “y - ‘SPh.r 

The definition (2.23) for [pjjj then allows (2.24) to be expressed in the 
form [p^J = [T^J[u Y J. 

The algebraic equations which result from the balance of traction force 


conditions between cell sides have been studied in [1]. Standard iterative 
methods appear to work quite successfully. 
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ii) A condensation method: 

The representation (2.8) for the transmission matrix for a triangular cell 

relates the traction values [p^J on the sides of a triangular cell 0 ) to 

the corresponding displacement values [v^J on the sides. As a result, with 

reference to Figure 6, j + j >P^_j ^ j »P^ + j j»P^ j_i can expressed in terms 

of the values v ±> j+ i »Vl , j ’ v i, j-1 * v i+l , j as wel1 as Using (2.13), 

we can also solve for v in terms of v and v . The result of 

1 »J i,J±l 

this condensation process, the details of which we shall not present here, is 
again the compact scheme (2.24). 


More generally , this process leads to a compact scheme for any cell based 
upon a condensation procedure for triangular cells, i.e., for the 
nonconforming approximation (2.5) condensation preserves compactness . 

An explanation of this is! Recall that for a triangle, the components of 
[v^J in (2.7) are the values 


[v y J = [v(_5 1 )Ay 1 ,vU 2 )Ay 2 ,vU 3 )Ay 3 J. 


(2.25) 


Because a linear nonconforming element is employed, each component is the 
value of the integral 


/ vdY ± . 
T i 


(2.26) 


In contrast, when using the conforming approximation (2.19) a trapezoidal 
approximation to (2.26) using vertex values is natural. In the latter case, 
however, the imposition of the balance of "traction' face condition does not 
reduce the number of displacement unknowns; this is in contrast to the 
nonconforming approximation (2.5) and highlights a basic difference between 
these two approaches. 
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II. 2. Compact Schemes for Elastic Bodies 

In write e and t in vector form as 


£- ( e ii> e 22 ,e 33 ,2e 23’ 2e 13 ,2e 12^ ’ 


1” ( t 11» t 22’ T 33’ 2t 23’ 2t 13 ,2t 12^ * 


(2.27) 


Let 3 = | — and 3 = (3,,3~,3_). Define symmetric matrix operators 

i 3x^ — 1 Z 3 

D l(3_),D 2 (8_) by 


3 0 0 


D 1 (3_) = 0 3 2 ° I ; 


0 0 3, 


0 9 o 

3 2 


°2 ( - " I 3 3 ° 3 i I » (2 * 28) 


3 3 0 

2 1 


and let 


D(3) = 


D lO) 


D 2 (3) 


with u = (u^u^u.j) 1 '. Then (1.1) can be written 


e(u) = D(3)u 


t(u) = Ce(u) 


where C is a symmetric 6x6 positive definite matrix. We shall restrict our 
discussion to the equilibrium problem 
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div t(u) = D C (a)x(u) = 0 


U = 


-r 


in ft, 
on r, 


(2.31) 


(i.e., q = 0 in (2.4)). 

Let n(£) = (n^ (O »i» 2 (£) >n^(0 ) denote the outward unit normal at a point 
L which lies at the center of a face y(J[) on the boundary y (oj) of a 

ce ll* If i* is a solution of (2.31) the traction force p(£) can be 

calculated as 


2.(0 = D t (n(C))*CDO)u 


(2.32) 


On a cell a) having k faces whose centerpoints are ***>?(> 

suppose are a basis set of solutions of (2.31) 

satisfying = 6 . Then 

k 

w(x) = l $>(x,E. )v(E. ) (2.33) 

i=l - - 1 - 1 

provides a solution of (2.31) which interpolates the constraint values v( £ ) 

— —l 

on y(o)). Recalling the notation 


let 


[ v y ] = [i(I 1 > A Y 1 ,vU 2 )Ay 2 , •••,v(C 4 )Ay £ ] 


(2.34) 


^(x;^) = *(x;5 1 )/Ay 1 . 


(2.35) 


Then, from (2.32) 




- | (D t (n(5))CD(3)>F( x ,5 )) .v(5,)4y 

i=i 1 *rl — 1 


(2.36) 


so that if 
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T Y<£»Ii> H Dt (n(5))CD(3)'F(x, il )| x!a5 

T y (£) = (2.37) 

then 

[tIJJ = (2.38) 


is the transmission matrix and 


P y (£) " T y (£)[v y J. (2.39) 

When k = 4, i.e., the case of a tetrahedron, the linear basis (l,x^,X2,x^) 
allows a simple construction. When k 2 4, the practical development of a 
basis set of equilibrium solutions in (2.33) can present a formidable problem 
when the constitutive matrix C in (2.30) is general. For k = 6 the 
result given in [2] can be applied directly when the cell is a brick; we shall 
describe this below. Otherwise the condensation technique indicated earlier 
can be used. Hence we can confine ourselves to examining the development for 
a tetrahedron. 

The Tetrahedral Cells 

Let (£^) indicate the centerpoints of the faces (y^) of a tetrahedral 
cell a); Aw is the volume and (Ay^) the areas of its faces. Similar to the 
example given earlier, use the centroid of the cell as origin and introduce 
volume coordinates for the tetrahedron (j^ ,5^ >^3 where 
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*T + AY i ln t a i ).x)/Aa 5 . 


(2.40) 


Then 


w(x) = 1 £(x,.L )v(£. > 

i=l 


(2.41) 


provides for the construction (2.33). 
Define (compare to (2.8)) 


T*(x,£) = (Aw) 1 D t (n(x))CD(n(£)), 


(2.42) 


Using (2.36)-(2.39) we obtain 


^ s T X )-[v r J - 


(2.43) 


This provides a compact scheme on w. 

Consider two tetrahedra w.aj' having the face Y ^ in common. The 

balance of traction condition at £ is 


pi(ii> + pi-<ii> = o. 


(2.44) 


This leads, using (2.43), to 


h h 4 h 

lvll.il> + y' (— ! 1 ) 1 ) - - l 2 T Y<Il»ij>I<Ij> A ^ 


“ J 2 V ( -i ’-P- ( -j )AY j 
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which may be solved for \r(£^) in the form 




j 2 + v«i- 


5pv(i')4rjJ. 


(2.45) 


2 

The equations should be compared to a similar result for V u = 0 given 
by (2.13). We believe, but have not been able to show, that the positive 
definiteness of the constitutive matrix C as well as geometrical properties 
arising from u) and induce properties in the coefficient matrices 

S y ,S^ in (2.45) which can allow standard iterative methods to be applied 
effectively. 


Brick Cells 

In order to develop a compact scheme for a brick cell, we could employ a 
condensation technique using the result just given for a tetrahedron. 
However, as indicated earlier, in this case a compact scheme has been derived 
from a slightly different viewpoint in Phillips and Rose [2]. Here we will be 
satisfied to present their scheme in a form closer to the notation presented 
earlier in this paper. 

Similar to the finite difference notations used in (2.24), let 


6 1 9 i,j,k “ 
P 1 = 


(0 i+V 2 ,j,k 9 i- V 2 ,j ,k )/2h i 

(9 i+V 2 J,k + 9 i-V 2 J,k )/2 


using similar definitions for >1*3 >63* Let ^ aces of w be defined 

by the planes x^ = (i ± V 2 )Axj , X2 = (j ±V 2 )Ax2> x^ = (k iV^Ax^, where 
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Ax v = 2h y ,v = 1,2,3. If n^ = (l,0,0) t ,n 2 = (0,1, 0) t , = (0,0, l) t , then the 

outward normal vectors on the faces are ±n v , v = 1,2,3. We let £^(±n y ) 
denote the traction force on the face for which in^ is the outward normal. 
Finally, again write 6_ = ^1*^2 ^3^* 

The compact scheme given in [2] can, with a suitable choice of certain 
nonessential parameters, be written as 

p h (±n) = ± D fc (n )CD(3)u + hAu u - X) v = 1,2,3 (2.46) 

— — v — v — — v v — — 


where 

a - 1 1 C % “)/( i O 

V=1 V=1 

using the notation indicated in (2.28)-(2.29) . 

The balance of flux condition across the face x v = const, common to 
neighboring cells u),u> is, then, simply 

° = fj = An^) + p h (-n y ) = D t (n v )CD(^)u - D t (n_ v )CDO)u 

+ h v 1 (p y (u + u ) - (X + X)) (2.47) 

In [2] a plane stress calculation involving an isotropic material was 
performed satisfactorily using several standard iterative methods. We refer 
the reader to that paper for details. 
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II. 3. Homogenization 

In the treatment of bonded materials, it is often useful to replace the 
laminate by an approximately equivalent homogeneous material. The methods 
considered in this paper offer an algebraic approach to this problem which we 
shall briefly describe. 

With reference to Figure 7, consider two bonded thin materials L and 
L' within which to and to' are typical small rectangular cells. We wish to 
treat to and to' as a single cell to with its own characteristic 

IV 

properties. Specifically, we wish to describe a compact scheme for to on the 
basis of the material properties reflected in the transmission matrices for 
to and to'. In order to accomplish this we exploit the principle that 
condensation preserves compactness. 


L 

h 

i 

OJ 


L' 

h' 


of 




6 



6 



2' 

• L 0 1 

'4 

a' 


b 


r 

• t o' < 

'3 

0) 



5 



W 

5 



Figure 7. Homogenization of a laminate. 
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Ass urae transmission matrices T and T , are known. As described 

to (0 

earlier, we may use a condensation technique to express the tractions on those 
sides of <o and to' which lie on the boundary of to in terms of the 
displacements which correspond to these sides. With reference to the figure, 
set Vj = v 2 = v_J 2, Vg = = v^/2 and, in terms of v^, and _v,_ 

v^, calculate p^ = (p^ + )/2 and _p^ = + £^)/2. We thus obtain a 

compact scheme for to which relates the tractions Op »P. »Pc»P A ) to the 

3 D “j D 

p values (v^ *^b *— 5 ’^6 ) and for which the transmission matrix 

to 

incorporates the material properties of both L and L ' . By this means, we 
suggest, L and L' can be treated as if it were a single material L. 
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